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1 1 ABSTRACT 

12 Conventional gravity and magnetic anomaly continuation invokes the standard Poisson boundary 

13 condition of a zero anomaly at an infinite vertical distance from the observation surface. This 

14 simple continuation is limited, however, where multiple altitude slices of the anomaly field have 

15 been observed. Increasingly, areas are becoming available constrained by multiple boundary 

16 conditions from surface, airborne, and satellite surveys. This paper describes the implementation 

17 of continuation with multi-altitude boundary conditions in Cartesian and spherical coordinates 

18 and investigates the advantages and limitations of these applications. Continuations by EPS 

19 (Equivalent Point Source) inversion and the FT (Fourier Transform), as well as by SCHA 

20 (Spherical Cap Harmonic Analysis) are considered. These methods were selected because they 

21 are especially well suited for analyzing multi-altitude data over finite patches of the earth such as 

22 covered by the ADMAP database. In general, continuations constrained by multi -altitude data 

23 surfaces are invariably superior to those constrained by a single altitude data surface due to 

24 anomaly measurement errors and the non-uniqueness of continuation. 

25 
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31 


I. INTRODUCTION 


32 In practice, no gravity or magnetic observation can be made with infinite accuracy. Thus, the 

33 further apart the truncated observations are from each other, the more biased they will be to 

34 different spatial attributes of the sources. This result makes it extremely problematic to relate, 

35 for example, truncated satellite altitude anomaly observations over hundreds of kilometers to 

36 imperfectly measured observations at the earth’s surface by simple downward continuation of the 

37 satellite data or upward continuation of the near-surface data (e.g., Schnetzer, 1985; Grauch, 

38 1993; von Frese et ah, 1999; 2005; Ravat et al., 2002; Kim et al., 2004; 2005). Measurement 

39 errors alone restrict the utility of spacebome anomaly continuations typically to within 

40 approximately ±50 km of mission altitude (von Frese et al., 2005). The problematic attributes of 

41 downward continuation, in particular, require that analyses of spaceborne geopotential 

42 observations be conducted essentially at satellite altitudes, unless they are integrated into a model 

43 that is also constrained by lower altitude observations. 

44 Measurement errors also result in uncertainties in the spectral properties of the anomaly 

45 observations that grow with increasing difference in the survey altitudes. Figure 1, for example, 

46 compares the spectral properties of North American crustal magnetic anomalies mapped by 

47 typical satellite surveys at 400 km (e.g., Magsat and CF1AMP) and the Decade of North 

48 American Geology (DNAG) compilation of aeromagnetic survey data (Hildenbrand et al., 1996). 

49 The comparison reveals a significant spectral gap between the two data sets with wavelengths 

50 comparable in scale to the magnetic effects of major crustal features. Due to measurement errors 

51 and the fundamental non-uniqueness of potential field continuation, there ultimately is no 

52 substitute to actual surveying for recovering these intermediate wavelength anomalies. Here, the 

53 intermediate wavelength anomalies might be substantially recovered by additional surveying at 

54 20 km altitude with U2 surveillance or other high-altitude aircraft (Hildenbrand et al., 1996). 

55 Global spherical harmonic models of the Earth’s gravity (e.g., Lemoine et al., 1998) and 

56 magnetic (Maus et al., 2009) fields are becoming increasingly available which are constrained by 

57 both satellite altitude and near-surface (i.e., terrestrial, marine and airborne) measurements. 

58 However, in contrast to the nearly global coverage that the satellite surveys provide, the coverage 

59 from near-surface surveys is far from global. Thus, when implementing near-surface predictions 


2 



60 from these models, particular care must be taken to verify that the predictions are constrained by 

61 actual survey data. 

62 Most geological gravity and magnetic anomaly studies tend to invoke local patches of 

63 coverage and complementary subsurface geological and geophysical constraints that are difficult 

64 to accommodate at the global scale of spherical harmonic modeling ( e.g ., von Frese and Kim, 

65 2003). Alternate approaches that are well suited for local applications in Cartesian and spherical 

66 coordinates include equivalent point source (EPS) inversion, the Fourier transform (FT) in 

67 Cartesian coordinates, and spherical cap harmonic analysis (SCHA) in spherical coordinates. 

68 These procedures are all applicable to gravity and magnetic studies of the Antarctic lithosphere 

69 south of 60°S. The sections below describe their uses for continuing multi-altitude data and 

70 consider the relative advantages and limitations of these applications. 

71 

72 II. EQUIVALENT POINT SOURCE (EPS) CONTINUATION 

73 In this section, EPS inversion is investigated for modeling anomaly continuations based on 

74 single- and dual-altitude boundary conditions and using satellite crustal anomalies to fill in the 

75 gaps in the near-surface survey coverage of the ADMAP database. EPS inversion has a long 

76 history of estimating interpolations, continuations, and other derivative and integral components 

77 of anomaly observations (e.g., Dampney, 1969; Ku, 1977; Mayhew, 1982; von Frese et al., 1981; 

78 1998; Asgharzadeh et ah, 2007; 2008). 

79 The EPS effect is the fundamental integrand for all forward modeling integrals in gravity and 

80 magnetic analysis (e.g., Blakely, 1995). For typical crustal ADMAP applications, it is 

81 convenient to express the gravity point pole and magnetic point dipole effects in terms of the 

82 spherical prism. More specifically, the gravity effect of the crustal prism in spherical (r, 0, cp) 

83 coordinates with uniform density (= mass/unit volume) contrast /I cr (e.g., von Frese et ah, 1981; 

84 Asgharzadeh et ah, 2007) is 

«/ >i ni d \ PR 

85 Ag(r, 0, (p) « A(p\ — Af7]4r. 2 )^.sin^.}4 , (1) 

i=\ /= i /= i oK K or 
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86 where G is the universal gravitational constant, R is the displacement between the source and 

87 observation points with respectively primed and unprimed coordinates, and (A„ Aj, A/) are the 

88 Gaussian-Legendre quadrature coefficients for the (nlxnjxni) point sources used to approximate 

89 the prism’s volume (e.g., Stroud and Secrest, 1966). In addition, 

90 Af>; -[(^ [(f -%)/2],A,; =[(r; -4,)/ 2] . (2) 

91 where ((p'i a , (p'ib), ( 9' ja , &jb), and (r' ia , r' lb ) are the lower (a) and upper (b) boundaries of the prism 

92 in the /-th coordinate of longitude (cp), the y'-th coordinate of co-latitude (6), and /- th radial 

93 coordinate (r). For completeness, note that Ku (1977) gives the equivalent Gauss-Legendre 

94 quadrature (GLQ) gravity effect of the prism in Cartesian coordinates. 

95 The total magnetic effect ( AT) of the crustal prism, on the other hand, with uniform 

96 magnetization (= dipole moment/unit volume) contrast Am in spherical coordinates (e.g., von 

97 Frese et ah, 1981; 1998; Asgharzadeh et ah, 2008) is 

nl nj ni i 

98 A7-(r,f?.p)*Ap',;E>f sin ^ , (3) 

/=1 7=1 1=1 A 

99 where the unit vectors £7' in source coordinates (r $’ <p’) and u in observation coordinates (r, 6, 

100 (p) are dotted into the source and observation coordinate gradient operators V' and V, 

101 respectively. The magnetization contrast Am represents the integrated effect of remanent and 

102 induced components, where the latter is the prism’s volume magnetic susceptibility times the 

103 intensity of the applied magnetic field at the source point. Again for completeness, note that Ku 

104 (1977) gives the equivalent magnetic effect of the prism by GLQ integration in Cartesian 

105 coordinates. 

106 Now, continuation, like any potential field analysis, can be generalized in matrix notation by 

107 AX = B, (4) 

108 where B is the (nxl) column matrix containing the n-anomaly observations, and AX is the 

109 forward model involving the (n x m) known coefficients of the design matrix A and the (mxl) 
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110 column matrix containing the m-unknown coefficients of the solution matrix X. The objective 

111 for continuation is to find the coefficients of X such that AX closely models B whereupon AX 

1 12 may be used evaluate the attributes of B at other coordinates. Most typically, the least squares 

1 1 3 solution is found by 

114 X = [A t A] _ 1 A t B (5) 

115 so that the sum of the squared residuals between the observations and the forward model’s 

116 predictions is minimum (e.g., von Frese et ah, 1988). 

1 17 Eq. 2 and Eq. 3 provide the forward models for relating respective gravity and magnetic 

1 1 8 anomalies to prism models that can be used for anomaly continuation. In particular, by 

1 1 9 specifying the geometric attributes of the prism model, the inversion (Eq. 4) can be set up to 

120 solve for the respective density and magnetization contrasts (Eq. 5) so that the modeled 

121 predictions (AX) match the observations (B) in a least squares sense. To obtain the 

122 coefficients of the design matrix A, the model is initialized by simply running it with the 

123 physical property (i.e., A a or Am) set to unity. 

124 In practice, however, errors in computing the A-coefficients due to the computer’s limited 

125 working precision and other uncertainties in the forward modeling may yield an unstable solution 

126 X with large and erratic values that are useless for predicting anything other than the original 

127 observations in B. To obtain a more stable and better performing solution, the system in Eq. 3 is 

128 commonly evaluated for the damped least squares solution 

129 X = [A t A + (EV) xl] _1 A T B , (5) 

130 where I is the identity matrix, and the scalar EV is variously called the damping factor, 

131 Marquardt parameter, or error variance (e.g., von Frese et al., 1988). The damped least squares 

132 approach requires choosing a value of EV that is just large enough to stabilize the solution for 

133 meaningful objective predictions (e.g., anomaly continuation, interpolation, etc.), yet small 

134 enough that the initial predictions provide an acceptable match to the observations B. For any set 

135 of observations and forward model, an essentially optimal value of EV can be selected from 

136 trade-off diagrams that contrast these two sets of model predictions for varying EV-values (e.g., 

137 von Frese et al., 1988). 
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138 Adapting the above results for anomaly continuation constrained by one or more boundary 

139 conditions is straightforward. Relating both airborne and satellite anomaly observations to a 

140 crustal prism model, for example, involves constructing the design matrix A = [A ajr bome A sate iiite] 

141 for the observation vector B = [B a jrbome B sa teiiite] T , where the submatrices Airborne and A sat eiiite 

142 reflect the geometric relationships between the crustal prism source coordinates and the 

143 coordinates of the airborne and satellite observations in the respective B a j r b 0 me and B sate iiite 

144 sub vectors. 

145 To test the performance of single- versus multi-altitude boundary conditions in anomaly 

146 continuation, Eq. 2 and Eq. 3 were used to model the respective gravity and magnetic effects of 

147 the five crustal prisms (Figure 2) outlined in top left anomaly maps of Figure 3 and Figure 4, 

"3 

148 respectively. The prisms were modeled for density contrasts ranging between -1.2 g/cm and 1 .8 

149 g/cm 3 , and cgs-magnetic susceptibility of contrasts of 0.0061468 and -0.028831 relative to the 

150 surrounding crustal rocks. All prisms were 25 km thick with tops at 5 km below sea level and 

151 located in Balkan region between 34° - 43° N and 21°- 30° W. the panel A of maps in Figure 3 

152 shows the modeled gravity effects to the nearest mGal over a 35^35 grid spanning the study area 

153 at altitudes of 5 km (top), 10 km, 100 km, 200 km, 400 km, and 600 km (bottom). Panel A of 

1 54 Figure 4 gives the complementary differentially reduced-to-pole (DRTP) magnetic effects to the 

155 nearest nT (e.g., von Frese et ah, 1981). The DRTP effects were obtained from Eq. 3 assuming 

156 vertical inclination of the applied field at all source and observation points with the field 

157 intensities at all source points taken from the World Magnetic Model (WMM) 2005 (NOAA, 

158 2005). 

159 The continuation analysis focused on relating the simulated airborne and satellite altitude 

160 effects at ,10 km and 400 km above sea level, respectively, to 50x50 arrays of crustal point 

161 masses and magnetic dipoles at depths of 40 km and 20 km, respectively, below sea level. These 

162 EPS inversions used the simpler forward model given by the integrands within the square 

163 brackets of Eq. 2 and Eq. 3 to obtain the relevant point mass (= Act x unit volume) and dipole 

164 moment (= dmxunit volume) solutions, respectively. In the two figures, the B-panels show the 

165 predictions from the EPS solutions constrained by the anomalies modeled at both 10 km and 400 

1 66 km above sea level, whereas the C- and D-panels give the continuations derived only from the 

167 single-surface effects at the respective 10 km and 400 km altitudes. Comparison of the dual- 
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168 surface continuations in the B-panels with the equivalent single-surface estimates in the C- and 

169 D-panels amply demonstrates the advantages of multi-altitude constrained continuation models 

1 70 especially for downward continuation. 

1 7 1 The comparisons at 5 km, and 1 0 km through 600 km altitudes at 1 0 km intervals are detailed 

172 further in Figure 5. The respective left and right panels demonstrate the performances of the 

173 gravity and magnetic continuations in terms of their correlation coefficients and root-mean- 

174 squared (RMS) differences with the effects of the five prisms. These results, of course, like those 

175 from any potential field analysis, are not unique. However, they clearly demonstrate the 

1 76 limitations of downward continuing satellite-only potential field observations. They also suggest 

177 that the multi-altitude magnetic continuations could be improved by incorporating additional data 

178 from around 60 km to 80 km altitudes. 

179 EPS models constrained by satellite and near-surface magnetic data are effective for filling in 

180 regional gaps in the near-surface survey coverage in the ADMAP compilation, which are 

181 particularly extensive in East Antarctica. The first compilation used satellite magnetic 

182 observations from the 400-km altitude, 6-month Magsat mission (Golynsky et al., 2001). 

183 However, these satellite data were collected during austral summer and fall periods, and thus 

184 maximally corrupted by external magnetic field variations. The 0rsted and CHAMP missions 

185 obtained much cleaner Antarctic satellite magnetic observations from several austral winters at 

186 altitudes of about 600 and 400 km, respectively, with measurement accuracies exceeding 

187 Magsat’s by roughly an order of magnitude. The initial Magsat-based predictions in the 

188 ADMAP compilation have been replaced by 0rsted- and CHAMP-based gap predictions (von 

1 89 Frese et ah, 2008). 

190 To obtain effective estimates of near-surface anomalies in the coverage gap, the inversion of 

191 the satellite data was optimized to match the near-surface anomaly observations around the gap’s 

192 perimeter (Kim et ah, 2004; 2007). Figure 6 illustrates the relative utility of this approach for a 

193 simulated regional gap in ADMAP’s near-surface survey coverage of the West Antarctic 

194 Peninsula and surrounding marine areas (Kim et al., 2004). Map A shows the original ADMAP 

195 10-km grid at 5-km above sea level of marine and airborne survey data for the well covered 

196 region. Due to the mission altitudes, the satellite data are mostly sensitive to the 400-km and 
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197 larger wavelength components of the near-surface data in Map B. The low-pass filtered data 

198 were related to a crustal prism model with negligible error that was evaluated at 400 km altitude 

199 with an accuracy of 3 nT in Map C to simulate the effects of Magsat measurement errors. The 

200 crustal prism model was also evaluated to a 0.3 nT accuracy at altitudes of 400 km in Map E and 

201 700 km in Map G to simulate the respective effects of CHAMP and 0rsted satellite observation 

202 errors. 


203 The simulated Magsat, CHAMP, and 0rsted data were combined with the near-surface data 

204 around the white-outlined periphery of the simulated gap in Map B and related to crustal prism 

205 models by least squares inversion. The resultant prism models were then evaluated for the 

206 respective Magsat-, CHAMP-, and 0rsted-based gap predictions shown in Maps D, F, and H. 

207 Table 1 compares the gap predictions against the original gap values from Map B in terms of the 

208 RMS differences, correlation coefficients (CC) and noise levels from the noise (N)-to-(S) signal 

209 ratio approximation (Foster and Guinzy, 1 967) given by 


210 


S /ICC I 


( 6 ) 


21 1 The comparison also includes gap predictions from the application of minimum curvature (MC) 

212 to the near-surface survey data. In general, the simulated results suggest that the use of the 

213 0rsted and CHAMP data improves the gap estimates by nearly 75% relative to the Magsat-based 

214 predictions. 


215 An opportunity to test the multi-altitude gap prediction method resulted when a new 

216 aeromagnetic survey was obtained in the regional gap located east of Coats Land and the 

217 Shackleton Range as shown in Map A of Figure 7 (Kim et ah, 2007). Map B gives the regional 

218 200 km and longer wavelength components of the new survey data that are much more consistent 

219 with the CHAMP- and 0rsted-based gap predictions in Map C than with the minimum curvature 

220 estimates in Map D and the near-surface predictions of the spherical harmonic geomagnetic field 

221 model MF4 (Maus and Rother, 2005) in Map E. 


222 Of course, as for any inversion, the results of continuation are not unique. Thus, they do not 

223 obviate the need for anomaly mapping in regions lacking measurements to better define the 

224 geologic implications in existing anomaly observations. 
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226 III. FOURIER TRANSFORM (FT) CONTINUATION 
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Gravity and magnetic anomalies are mapped at elevations accessible by borehole, tunnel, ground 
and ice surface, ocean and lake bottom, submarine, ship, airplane, balloon, space shuttle tether, 
and satellite surveys. Where the horizontal dimensions of the survey are several hundred 
kilometers or less, the data are commonly gridded at constant elevation and horizontal intervals 
measured in linear units of length for analysis. In this case where the anomaly data are registered 
to a Cartesian grid, effective anomaly continuation can be done most elegantly and efficiently 
using the Fourier transform as the forward model of the inversion (e.g., Blakely, 1995). 

The free-space behavior of gravity and magnetic anomalies S is governed by Laplace’s 
equation 


/ a 2 d ! 

( r b 

dx 2 dy 2 



= 0 


(7) 


so that the vertical z-dimensional properties of the anomalies may be assessed from their mapped 
x- and/or y-dimensional attributes. In particular, for gridded S(x, y) with the gridded Fourier 
transform S(I,k) in wavenumbers / and k, taking the Fourier transform of Eq. 7 gives 


d 2 s 


dz 2 


= (2/r)T//+ f t 2 )S > 


( 8 ) 


where (fi = l/M) and (fk = k/N) are the linear frequencies for the number of anomaly values M and 
Ain the x- and y-directions, respectively. Eq. 8 is a second order differential equation with 
constant coefficients that has the general solution 

S(z)=Axe +:fa +Bxe :fa , (9) 


where 

a = 4;r 2 (/ / 2 +//) (10) 

Boundary conditions are commonly invoked to specify the form of the solution in Eq. 9 that 
is appropriate for application. By the classical Poisson boundary condition, acceptable solutions 
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are those that approach zero as the elevation z becomes infinite. Thus, for upward continuation 
where z > 0, the boundary condition requires A = 0 so that the solution becomes 

= ( 11 ) 

Referencing the grid to z = 0 shows that B - S(0) = S so that the upward continuation of 
the anomaly grid to any positive z > 0 is 

S(z\ paN -Sx e -=' r \ (12) 

Similarly, for downward continuation where z < 0, the boundary condition requires B = 0 so 
that the solution is 


Stf^Axe ^ . (13) 

Again, taking the elevation of the grid as relative zero shows that A = S(0) = S so that the 
downward continuation of the anomaly grid to negative z < 0 is 

SM^Sxe^’. (14) 

Inversely transforming Eq. 12 and Eq. 14 is the conventional approach for upward and 
downward continuing a single anomaly grid through positive and negative elevations about its 
elevation (e.g., Pawlowski, 1995; Fedi and Florio, 2002; Cooper, 2004). The singular uses of the 
A and B coefficients in Eq. 9 define the two end member solutions for downward and upward 
continuation, respectively. However, multiple altitude anomaly grids are becoming increasingly 
available where both coefficients may be applied to evaluate the anomaly values at the 
intermediate altitudes between the grids. 

Specifically, let the lower altitude anomaly grid be referenced to relative zero and the higher 
altitude grid be referenced to an elevation H where the Fourier transforms of S(0) and S(H) are 
5(0) and S(H ) , respectively. Both transforms may be related through Eq. 9 to common A and 
B coefficients given by 


md B = s i H 1 

e' zr ° 


-S( 0)e +H ^ 




(13) 
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280 


so that the continued field at an intermediate z-altitude for all 0 < z < H is 


281 

282 = i_' 

t 

283 

284 

285 with 

286 

287 

288 Figure 8 illustrates the use of the dual Fourier transforms for two altitude slices of the noise- 

289 free gravity effects for a simple model of two prisms. The gravity effects of the prisms were 

290 calculated in Cartesian coordinates at relative 0 km (Map A) and at 1 6 km (Map B) and 20 km 

291 (Map C) above the relative 0-km effects. Map D shows the dual altitude continuation estimates 

292 at 16 km and Map E gives the difference between Map B and Map D. Map F is the 

293 conventional downward continuation of Map C at 16 km, whereas Map G is the 16-km upward 

294 continuation of Map A with Map H and Map I giving the related differences from the modeled 

295 effects in Map B. The dual and single surface continuations are largely consistent, although the 

296 dual surface estimates in Map D are clearly superior to the 4-km downward and 1 6-km upward 

297 continued single surface estimates in Maps F and G. 

298 The general consistency of the single- and dual-surface continuations in Figure 8 is mostly 

299 due to the source effects being minimally distorted in the two anomaly grids. This level of 

300 source signal fidelity is not likely to occur in practice where, for example, in the lower altitude 

301 grid the effects of near-surface sources mask within working precision elements of the more 

302 regional source effects that in the higher altitude grid are preferentially expressed at the expense 

303 of the shorter wavelength near-surface source effects. Measurement and anomaly reduction 

304 errors also will substantively affect the fidelity of gridded source signals at the different altitudes. 



11 



305 As an example. Figure 9 considers the altitude variations of the reduced-to-pole magnetic 

306 anomalies over the southern Greenland area constrained by airborne (Verhoef et ah, 1996) and 

307 CHAMP satellite (Maus et ah, 2002) measurements. The magnetic data were Fourier transformed 

308 assuming the flat-earth approximation given the approximate 1,600 km by 1,600 km area of the 

309 study (Strang van Hees, 1990). The single surface transform was used to upward continue the 

310 airborne data to 20 km in Map A, 50 km in Map B, 250 km in Map C, 300 km in Map D, and 

311 350 km in Map E. However, the CHAMP data at 350 km in Map F show a considerably 

312 stronger effect over Iceland than is apparent in the upward continued airborne data of Map E. 

313 Indeed, the dual surface FT continuations at 300 km in Map G, 250 km in Map H, and 50 km in 

314 Map I suggest that this disparity persists at altitudes of 250 km and higher. Clearly, the CHAMP 

315 data provide important additional constraints for developing a comprehensive crustal magnetic 

3 1 6 model of the study area. 

317 

318 IV. SPHERICAL CAP HARMONIC ANALYSIS (SCHA) 

319 Spacebome and lower altitude magnetic and gravity anomalies obviously provide different 

320 insights for resolving the crustal sources of the potential field anomalies. In particular, the lower 

321 altitude anomalies provide much more detail on the crustal sources than can be observed at 

322 satellite altitude. The satellite anomalies, on the other hand, typically yield better regional images 

323 of the sources than can be obtained from regional compilations of the lower altitude anomalies 

324 that have been mapped with disparate survey parameters. At regional scales of coverage such as 

325 represented by the ADMAP compilation for the Antarctic area south of 60° S, the spherical cap 

326 harmonic synthesis of satellite and near-surface potential field observations is particularly useful 

327 and elegant model for crustal modeling. For example, it is analytically more useful than the 

328 simple grid for representing the ADMAP compilation because it can directly estimate magnetic 

329 anomalies and their gradient and tensor components anywhere on and above the Antarctic 

330 surface (e.g., von Frese et ah, 2008). 

331 Spherical cap harmonic analysis (SCHA) was developed by Haines (1985) for geomagnetic 

332 main and external magnetic field modeling at regional scales. It has also been extended to 

333 regional magnetic studies of the lithosphere (e.g., Haines and Newitt, 1986; De Santis et ah, 

334 1990; Kotze, 2001; Thebault et ah, 2006; Thebault and Gaya-Pique, 2008; von Frese et al., 


12 



335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 


2008), as well as gravity studies (e.g., Hwang and Chen, 2007; De Santis and Torta, 1997). This 
approach adapts the global spherical harmonic solution to a spherical cap solution for the 
magnetic potential 


V{r,6,cp) 


' V im « 

II 

w=0 m~Q 


a(~ )' 

r 




( g"' J cos mcp+ h"' 1 s in in<p )P m n ( cos 0 ) + 


Y Y a(—) n (g m /cosmA,(p + ti"' e s\nm(p)p™(cosQ), 

n I nr () r 


(16) 


where N in , and N ex , are the maximum indices for internal and external sources, respectively, P„ m 
is a Legendre function of degree n and order m, and g™’ 1 , h™’* , g'”’ e , h™ ,e are the spherical 
cap Gauss coefficients. If the half-angle of the spherical cap is denoted byi9 0 , the rik(m) for a 
given m are determined as the roots of 

dp"!(cosd o )/d0 = O, k- m = even, (17) 

and additionally, where differentiability with respect to 6 is required, of 

dp™(cos6 0 ) = 0, k- m = odd. (18) 


Truncating the expansion in Equation (16) at k = K sets the number of model coefficients at 
(K+lf. When applying SCHA, a global spherical harmonic potential is typically removed from 
the total potential to improve convergence as well as extrapolation beyond the spherical cap 
boundary (Haines, 1985). In general, the maximum degree of a spherical cap harmonic model, 
truncated at a certain number of model coefficients, is a function of the half-angle of the 
spherical cap that determines the minimum wavelength resolution that can be obtained (e.g., 
Kotze, 2001). As a radial distance is also included in the formulation, data over a range of 
altitudes can be used to derive a model. 
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361 Two examples are presented in this study. The SCHA is a parametric method with the 

362 coefficients of the harmonic representation of the gravity and magnetic potential being model’s 

363 parameters. This method has been proven to work with both gridded and irregularly distributed 

364 data, which in example are along track data. We used a variation of the SCHA, in which the 

365 spherical cap is mapped onto a hemisphere, thus ensuring the proper harmonic continuation of 

366 the gravity and magnetic field anomalies and to different altitudes between satellite and the Earth 

367 surface. 

368 

369 Figure 10 shows the use of SCHA to model the effects of the NE Asia gravity field 

370 developed by the Curtin University of Technology (CUT) in Western Australia (Kuhn and 

371 Featherstone, 2005). The gravity anomalies were synthesized on a 1° x 1° grid with the center of 

372 the spherical cap is 50 0 N and 130 0 E. In this example, the synthetic gravity anomaly data at the 

373 Earth’s surface (c) and 20 km altitude (a) were jointly modeled by the SCHA with negligible 

374 error to obtain gravity anomalies at 10 km altitude (b) which show the compiled gravity 

375 anomalies at 10 km altitude with a resolution of 200 km half-wavelength. The model errors 

376 account for 10% of the gravity anomaly magnitude at most and can be further reduced by 

377 increasing the spectral resolution of the model. Table 2 shows the statistics of synthetic gravity 

378 data at 20, 10, and 0 km and correlation coefficients shown in Figure 10. 

379 

380 Another example modeling the radial geomagnetic field component at satellite (a), U2 (b) 

381 altitudes and the Earth surface (c) is given in Figure 11. Here, the satellite altitude (685 km) 

382 component for December 1999 was globally evaluated from the three-axis magnetometer 

383 measurements of the KOrea Multi Purpose Satellite-1 (KOMPSAT-1) that provided the 

384 satellite’s attitude control (Kim et al., 2007) as shown Figure 11 (a). These data were used to 

385 estimate the radial components by the SCHA at satellite (685 km), U2 (20 km) altitudes and at 

386 the Earth surface (0 km). Table 3 summarizes statistics and correlation coefficients of radial 

387 component of KOMPSAT-1 measurements estimated at satellite (685 km), U2 (20 km), and 

388 Earth surface (0 km) as shown in Figure 11. This example illustrates a further use of SCHA in 

389 developing improved regional core field models for extracting crustal components in the survey 

390 data. 

391 
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392 The Antarctic Reference Model (ARM), for example, was developed for ADMAP 

393 applications to best represent the spatial and temporal properties of the Antarctic core field as 

394 recorded over the period 1960-2002 by an international network of 28 observatories and satellite 

395 magnetic data from six missions (Gaya-Pique et al., 2006). The SCHA predictions from ARM 

396 have greater sensitivity for the Antarctic spatial and temporal variations of the core field than do 

397 conventional global spherical harmonic geomagnetic field models like the IGRF-10 (IAGA, 

398 2005). 

399 Data from the Arctowski, Scott Base, Syowa, and Vostok geomagnetic observatory data 

400 through 2003, for example, showed that the ARM estimates of the annual variations of the 

401 Antarctic core field are roughly 31% and 22% better than the IGRF-10 and CM4 (Sabaka et al., 

402 2004) estimates, respectively. In addition, the core field predictions were checked at the Georg 

403 von Neumayer observatory where data were obtained that had not been included in the original 

404 ARM model. At this observatory, which was operated over the period 1983.5-1991.5 on the 

405 coast of Dronning Maud Land (70.617°S, 351.633°E), the ARM estimates were found to be 76% 

406 and 75% better than the respective IGRF-10 and CM4 estimates. Additionally, ARM estimated 

407 the mean secular variation some 78% and 75% better than the respective IGRF-9 and CM4 

408 models. 

409 

410 V. CONCLUSIONS 

411 Low-earth-orbiting satellites are recording state-of-the-art magnetic and gravity field 

412 observations of the earth at altitudes ranging over roughly 300-700 km. However, the 

413 measurement errors and non-uniqueness of the continuation process severely limit estimating 

414 these anomalies at the near-surface. Indeed, numerical anomaly simulations from satellite to 

415 airborne altitudes suggest that effective downward continuations of the satellite data are 

416 restricted to within approximately 50 km of the observation altitudes while upward continuations 

417 can be effective over a somewhat larger altitude range. Given current anomaly error levels, the 

418 multi-altitude inversion of satellite and near-surface anomalies is a particularly promising 

419 approach for implementing satellite geomagnetic observations for crustal studies. 
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420 This study investigated multi-altitude implementations of equivalent point source (EPS) and 

421 spherical cap harmonic analysis (SCHA) modeling in spherical coordinates and the Fourier 

422 transform (FT) for continuation in Cartesian coordinates. These implementations revealed that 

423 the anomaly estimates at the intermediate altitudes from multi-altitude constrained models were 

424 invariably superior to the classical downward and upward continuations of a single slice of the 

425 anomaly field. In addition, near-surface downward continuations of satellite altitude 

426 observations or satellite altitude upward continuations of near-surface data can be severely 

427 problematic due to measurement and data reduction errors and the inherent non-uniqueness of 

428 any continuation. As a result, an observation’s capacity to define a field’s spatial behavior is 

429 restricted in practice to a relatively small region about the observation site beyond which 

430 mapping is required to confidently resolve the more distant spatial details of the field. In other 

43 1 words, there ultimately is no substitute for mapping to establish the value of an anomaly field at 

432 an unsurveyed location. 

433 
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561 gap in aeromagnetic anomaly coverage. The prediction statistics include the root-mean- 

562 square (RMS) difference, the correlation coefficient (CC), and related noise (N) levels. In 
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570 measurements estimated at satellite (685 km), U2 (20 km), and Earth surface (0 km) as 
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575 Figure 1. Normalized amplitude spectra from geomagnetic surveys at typical satellite (400 km), 

576 U2 (20km), and airborne (1 km) altitudes (adapted from Hildenbrand et al., 1996). 

577 

578 Figure 2. Location of the 5 crustal prisms. All prisms were 25 km thick with tops at 5 km below 

579 sea level and located in Balkan region between 34° - 43° N and 2 1°- 30° W. 

580 

581 Figure 3. Gravity effects (mGal) from five prisms in the spherical coordinates modeled (A) 

582 directly from the GLQ integration, (B) by joint EPS inversion of the 10 km and 400 km 

583 fields, (C) by single EPS inversion with the 10 km field, and (D) by single EPS inversion 

584 with 400 km field. 

585 

586 Figure 4. Magnetic effects (nT) from five prisms in the spherical coordinates modeled (A) 

587 directly from the GLQ integration, (B) by joint EPS inversion of the 10-km and 400-km 
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fields, (C) by single EPS inversion with the 10-kni field, and (D) by single EPS inversion 
with 400-km field. 

Figure 5. Comparison between the modeled and estimated (A) gravity (mGal) and (B) magnetic 
fields (nT) in Figure 3 and Figure 4, respectively: (upper) correlation coefficient, 
(middle) root-mean-square (RMS) difference, and (lower) Log (RMS). 

Figure 6. Estimating near-surface magnetic effects to fill in regional gaps in the ADMAP 
airborne and marine anomaly survey coverage (Kim et ah, 2004). (A) ADMAP near- 
surface magnetic anomalies for the West Antarctic Peninsula and surrounding marine 
areas at 5 km altitude. (B) Map A low-pass filtered for 400-km and longer wavelengths 
with simulated gap outlined in white. (C) Simulated Magsat anomalies at 400 km altitude 
from a crustal prism model of Map B evaluated to 3-nT accuracy. (D) Near-surface 
magnetic anomaly estimates at 5 km altitude with the coverage gap filled in using the 
simulated Magsat in Map C. (E) Simulated CHAMP anomalies at 350 km altitude 
evaluated at 0.3-nT accuracy. (F) Near-surface magnetic anomaly estimates at 5 km 
altitude with the coverage gap filled in using the simulated CHAMP data in Map E. (G) 
Simulated 0rsted anomalies at 700 km altitude evaluated at 0.3-nT accuracy. (H) Near- 
surface magnetic anomaly estimates with the coverage gap filled in using the simulated 
0rsted data in Map G. 

Figure 7. Comparisons of magnetic anomaly gap predictions with a recently surveyed portion of 
the regional gap located east of Coats Land and the Shackleton Range in Antarctica (Kim 
et al., 2007). (A) Recently surveyed aeromagnetic anomalies within the red-bordered 
rectangular area and (B) their 200-km and longer low-pass filtered components compare 
well with (C) the gap predictions over the survey area obtained using joint CHAMP and 
0rsted data and the earlier aeromagnetic anomalies. The gap predictions from (D) 
minimum curvature and (E) the spherical harmonic geomagnetic field model MF4 (Maus 
and Rother, 2005) are also shown. The correlation coefficients between the regional 
survey components in Map B and the predictions in Map C, Map D, and Map E are 
0.45, -0.5, and -0.18, respectively. 
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619 

620 Figure 8. Gravity anomalies calculated at 0 km (A), 16 km (B), and 20 km (C) from prism 

621 models. Map B gives the differentially continued anomalies at 16 km from the multi- 

622 altitude FT using the anomalies from Map A and Map C. Map B gives the related 

623 anomaly differences (B - D) with the (minimum, maximum), mean, and standard 

624 deviation values listed along the bottom. Map F and Map G were calculated by 

625 conventional single-surface Fourier transform as the downward continuation of Map C 

626 and upward continuation of Map A, respectively. Their anomaly differences relative to 

627 the modeled effects at 16 km in Map B are given in Maps H and Map I, respectively. 

628 

629 Figure 9. (A) Aeromagnetic anomalies reduced-to-pole for Southern Greenland and surrounding 

630 marine areas at 20 km altitude (Verhoef et al., 1996) and conventionally upward 

631 continued to (B) 50 km, (C) 250 km, (D) 300 km, and (E) 350 km. (F) CHAMP based 

632 MF6 model estimates at 350 km (Kim et al., 2010) that were combined with the anomaly 

633 data in Map A for the dual-surface FT continuations at (G) 300 km, (H) 250 km, and (I) 

634 50 km. Selected correlation coefficients between the maps are CC(E,F) = 0.16, CC(D,G) 

635 = 0.61, CC(C,H) = 0.88, and CC(B,I) = 0.99. 

636 

637 Figure 10. Spherical cap harmonic analysis (SCHA) model of joint gravity anomalies at 20 

638 km (a) and 0 km (c) altitudes evaluated at the intermediate altitude of 10 km (b). 

639 

640 Figure 11. Spherical cap harmonic analysis (SCHA) model of radial geomagnetic field 

641 components evaluated from global KOMPSAT-1 three-axis magnetometer 

642 measurements during December 1999 at 685 km altitude (a) and SCHA model 

643 estimates at the 20 km altitude of a U2 reconnaissance survey (b) and at the Earth 

644 surface (d) (adapted from Kim et al., 2007). 
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Table 1. Performance statistics for using minimum curvature (MC) and Magsat (Figure 3(B)), 0rsted 
(Figure 3(D)), and CHAMP (Figure 3(F)) magnetic anomalies to fill a simulated gap in aeromagnetic 
anomaly coverage. The prediction statistics include the root-mean-square (RMS) difference, the 
correlation coefficient (CC), and related noise (N) levels. In the right four columns, the upper triangular 
portion gives the relative noise reductions for the gap predictions from the various constraints. 
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Table 2. Statistics of synthetic gravity data at 20, 10, and 0 km and correlation coefficients shown in 
Figure 9 (unit: mGal). 



Min 

Max 

Standard 


Correlation Coefficient 


Deviation 

at 20 km 

at 10 km 

at 0 km 

at 20 km 

-284.87 

250.47 

65.01 

- 

0.998 

0.991 

at 1 0 km 

- 335.83 

280.13 

69.60 

0.998 

- 

0.997 

at 0 km 

- 409.66 

319.27 

75.27 

0.991 

0.997 
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Table 3. Statistics and correlation coefficients of radial component of KOMPSAT-1 measurements 
estimated at satellite (685 km), U2 (20 km), and Earth surface (0 km) as shown in FigurelO (unit: nT). 
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Standard 


Correlation Coefficient 



Deviation 

at 685 km at 20 km 

at 0 km 

at 685 km 

- 55378.50 

- 13201.80 

8821.38 

- 

0.991 

0.990 

at 20 km 

- 70959.30 

-22931.40 

10743.58 

0.991 

- 

1.0 

at 0 km 

- 71515.10 

- 23270.30 

10813.16 

0.990 

1.0 
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Figure 1. Normalized amplitude spectra from geomagnetic surveys at typical satellite (400 km), U2 
(20km), and airborne (1 km) altitudes (adapted from Hildenbrand et al., 1996). 
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Figure 2. Location of the 5 crustal prisms. All prisms were 25 km thick with tops at 5 km below sea 
level and located in Balkan region between 34° - 43° N and 21°- 30° W. 
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Figure 3. Gravity effects (mGal) from five prisms in the spherical coordinates modeled (A) directly 
from the GLQ integration, (B) by joint EPS inversion of the 10 km and 400 km fields, (C) by single 
EPS inversion with the 10 km field, and (D) by single EPS inversion with 400 km field. 
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Figure 4. Magnetic effects (nT) from five prisms in the spherical coordinates modeled (A) directly 
from the GLQ integration, (B) by joint EPS inversion of the 10-km and 400-km fields, (C) by single 
EPS inversion with the 10-km field, and (D) by single EPS inversion with 400-km field. 
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Figure 5. Comparison between the modeled and estimated (A) gravity (mGal) and (B) magnetic fields 
(nT) in Figure 3 and Figure 4, respectively: (upper) correlation coefficient, (middle) root-mean-square 
(RMS) difference, and (lower) Log (RMS). 





Figure 6. Estimating near-surface magnetic effects to fill in regional gaps in the ADMAP airborne and 
marine anomaly survey coverage (Kim et al., 2004). (A) ADMAP near-surface magnetic anomalies for 




the West Antarctic Peninsula and surrounding marine areas at 5 km altitude. (B) Map A low-pass 
filtered for 400-km and longer wavelengths with simulated gap outlined in white. (C) Simulated 
Magsat anomalies at 400 km altitude from a crustal prism model of Map B evaluated to 3-nT accuracy. 
(D) Near-surface magnetic anomaly estimates at 5 km altitude with the coverage gap filled in using the 
simulated Magsat in Map C. (E) Simulated CHAMP anomalies at 350 km altitude evaluated at 0.3-nT 
accuracy. (F) Near-surface magnetic anomaly estimates at 5 km altitude with the coverage gap filled in 
using the simulated CHAMP data in Map E. (G) Simulated 0rsted anomalies at 700 km altitude 
evaluated at 0.3-nT accuracy. (H) Near-surface magnetic anomaly estimates with the coverage gap 
filled in using the simulated 0rsted data in Map G. 




Figure 7. Comparisons of magnetic anomaly gap predictions with a recently surveyed portion of the 
regional gap located east of Coats Land and the Shackleton Range in Antarctica (Kim et al., 2007). (A) 
Recently surveyed aeromagnetic anomalies within the red-bordered rectangular area and (B) their 200- 
km and longer low-pass filtered components compare well with (C) the gap predictions over the survey 
area obtained using joint CHAMP and 0rsted data and the earlier aeromagnetic anomalies. The gap 
predictions from (D) minimum curvature and (E) the spherical harmonic geomagnetic field model MF4 
(Maus and Rother, 2005) are also shown. The correlation coefficients between the regional survey 
components in Map B and the predictions in Map C, Map D, and Map E are 0.45, -0.5, and -0.18, 
respectively. 
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Figure 8. Gravity anomalies calculated at 0 km (A), 16 km (B), and 20 km (C) from prism models. 
Map B gives the differentially continued anomalies at 16 km from the multi-altitude FT using the 
anomalies from Map A and Map C. Map B gives the related anomaly differences (B - D) with the 
(minimum, maximum), mean, and standard deviation values listed along the bottom. Map F and 
Map G were calculated by conventional single-surface Fourier transform as the downward 
continuation of Map C and upward continuation of Map A, respectively. Their anomaly differences 
relative to the modeled effects at 16 km in Map B are given in Maps H and Map I, respectively. 


Figure 9. (A) Aeromagnetic anomalies reduced-to-pole for Southern Greenland and surrounding 
marine areas at 20 km altitude (Verhoef et al., 1996) and conventionally upward continued to (B) 50 
km, (C) 250 km, (D) 300 km, and (E) 350 km. (F) CHAMP based MF6 model estimates at 350 km 
(Kim et al., 2010) that were combined with the anomaly data in Map A for the dual-surface FT 
continuations at (G) 300 km, (H) 250 km, and (I) 50 km. Selected correlation coefficients between 
the maps are CC(E,F) = 0.16, CC(D,G) = 0.61, CC(C,H) = 0.88, and CC(B,I) = 0.99. 
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Figure 10. Spherical cap harmonic analysis (SCHA) model of joint gravity anomalies at 20 km 
(a) and 0 km (c) altitudes evaluated at the intermediate altitude of 10 km (b). 
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Figure 11. Spherical cap harmonic analysis (SCHA) model of radial geomagnetic field 
components evaluated from global KOMPSAT-1 three-axis magnetometer measurements 



during December 1999 at 685 km altitude (a) and SCHA model estimates at the 20 km altitude 
of a U2 reconnaissance survey (b) and at the Earth surface (d) (adapted from Kim et ah, 2007). 



